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1 . 0 Introduction 

The problem of interpolating or approximating data from 
scattered measurements arises in many areas of science and 
engineering. The importance of the problem has resulted in a 
large number of methods for solution of the problem, as has been 
noted in surveys by Schumaker [55], Barnhill [4], and more re- 
cently by Franke [19]. The problem can be described very easily. 
Given data (x^,y^,f^), i=l,...,N, construct a (smooth) function, 
F(x,y) such that F(x^,y.) = f ^ , i=l,...,N. No assumption is made 
regarding the spacing of the independent variable data. More 
generally, especially when the data are subject to errors, one 
may wish to relax the interpolation condition and approximate the 
surface. The problem has an obvious generalization to more 
independent variables. The existence of many methods for such a 
surface is due to the many sources of data and the great number 
of possible dispositions of data points. Some advice on how to 
proceed for given data is contained in Sabin [50]. It is clear 
that no one method is satisfactory for all cases. 

In the last few years a number of advances have been made. 
Because of space and time limitations, this paper will be con- 
fined to a few of those which I feel to be particularly important 
or interesting. The first area concerns the mathematical under- 
pinnings of the multiquadric method of Hardy [22] . This method 
has previously been noted to work well in a variety of cases (see 
Franke [19] and Kansa [29], for example), but until recently 
little was known in terms of a mathematical theory for the 
method. These developments will be reviewed in section 2. The 
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construction of surfaces with tension parameters or satisfying 
constraints seems to be desirable based on the importance of such 
ideas in the univariate case, and recent results are discussed in 
section 3. Finally, since the amount of data is sometimes far in 
excess of what is required to define a sufficiently accurate 
surface, the problem of surface approximation and selection of 
subsets on which to base interpolating (or approximating) sur- 
faces has become important, and this topic will be treated in 
section 4 . 

There have been other developments which are important, but 
will not be treated here. These include convergence of Shepard’s 
method (see Farwig [17]), interpolation on the sphere (see Lawson 
[32], Wahba [63], Renka [49], Ramaraj [48], and Nielson and 
Ramara.j [46]), and some new implementations of "finite element" 
methods and related things such as derivative estimation and new, 
higher order elements (see Alfeld [2-3], LeMehaute [33], [34], 
Sablonniere [51]). Certain methods based on statistical ideas 
(especially Kriging in geology and Optimum Interpolation in 
meteorology) continue to be the focus of much effort in those 
particular disciplines. Relevant recent references include 
Journel [28] and Thiebaux [59]. 

2.0 Multiquadrics and Related Methods 

The multiquadric (MQ) method was proposed by Hardy [22], and 
he has investigated the fitting properties of the method when 
applied to data from various sources in a series of papers that 
have appeared since that time [23-25]. The scheme is quite 
simple to implement and is reasonably efficient in terms of 
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computer resources provided the number of points is not large. A 

2 2 1/2 

basis function (a quadric) B.=(d. +r)^ is associated with 

0 3 

•f* Vv 

the data point. Here d. is the distance from the point (x,y) 

V 

to the data point (x.,y.) and r is a parameter in the method. 

Note that each of the basis functions is a radial function with 
respect to the data point, and that they are all translates of 
each other, which overcomes the usual problem of running afoul of 
the Haar theorem regarding interpolation in more than one 
variable. Now, a linear combination of the functions 



N 

F(x,y) = B^(x,y) 

j = l 

is required to interpolate the given data. This yields t,he 
system of equations, 



N 



j = l 




i=l,2, . . . ,N 



The existence of the interpolant is dependent on the nonsingu- 
larity of the matrix {B.(x.,y.)}. Other authors have also inves- 

± 1 . 

tigated the scheme, generally from an empirical point of view 
[19], [27], [57]. These authors have found that the method is 
quite adept at yielding good approximations, in many situations 
better than any other method. In addition to this method, the 
use of the reciprocal of the above basis function leads to the 
"reciprocal" MQ method. 

The MQ method is one of a class that I have previously 
called "global basis function methods", and which others have 
called "kernel" methods. In general the approximation takes the 
form 
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N M 

F(x,y) = Bj(x,y) + ^bj Pj(x,y) . 

j=l j=l 

where {p.} is a set of monomials of degree <m. The equations for 

the coefficients in such methods may be written in the form 
N M 

^Bj(x.,y^)aj + = 'i ■ i=l N . 

J=l 0=1 



N 

2^Pi(Xj,yj)aj = 0 . i=l M . 

J = l 

The first of the equations require interpolation to the data by a 

linear combination of the basis functions B^(x,y) plus the M 

polynomial terms, P-(x,y), while the last set requires the coef- 

3 

ficients to satisfy a certain constraint, which as we shall see, 
may related to the conditional positive definiteness of 
{B^(x^,y^)}, and serves to guarantee exactness for the set of 
polynomials {Pj(x,y)}. In matrix form, we may write 





2,1 The Multiquadric Method: Theory 

The intriguing aspect of the method is that until recently, 
very little was known in terms of a mathematical basis for the 
efficacy of the method, even whether or not the coefficient 
matrix was possibly singular. As recently as 1983 at The Inter- 
national Symposium on Surface Approximation, in Gargnano , Italy, 
I proposed as a conjecture the inequality 
(-1)^ ^det{Bj (x^ ,y^ ) } > 0 . 

As it turned out, Charles Micchelli promptly heard of the conjec- 
ture, and subsequently proved it, and along the way, theorems 
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which answered some other questions as well. 

To discuss the results of Micchelli [41] it is necessary to 
define some terms and give some background information. Because 
of the generalization to s-dimensional space, for this discus- 
sion, points in R® will be denoted as (possibly superscripted) 
boldface vectors rather than subscripted coordinates. 

Definition: A continuous function F(t), defined on [0,oo) is said 

to be conditionally (strictly) positive definite of order k on R 

if for any distinct points x"*" , x'^ , . . . , x” in R^ , and scalars c^ , 

Cg, ...» c^ such that 
n 

^^CiP(x^) = 0 

i = l 

for all polynomials p over R of degree <k, the quadratic form 
i=l j=l 

is (positive) nonnegative. 

Let the class of functions which are conditionally positive 

o C! 

definite of order k over R be denoted by Pj^(R ) > ^nd the class 
of functions which are conditionally positive definite of order k 
over R^ , for all s, by P^. Further, recall that a function F is 
completely monotonic on (0,co) if it is in C (0,co) and 
(-1 )"‘f("‘) (t) > 0 for all t>0 and m=0,l,2 

Theorem 1: F is in whenever F is continuous on [0,co) and 

(-1) F' '^(t) is completely monotonic on (0,oo). 

This theorem is due to Schoenberg [54] for order k=0 . Fur- 
ther, if F is not a polynomial of degree <k, and if the points in 
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the definition of completely monotonic are distinct, then the 
quadratic form is positive. Take F(t) = (r+t) ' , and note 
that for m> 1 , ( t ) =C^ ( -1 )'^ ^(r^+t)^^ 2m)/2^ where is a 
positive constant. This special case leads to the coefficient 
matrix for the MQ method, which is seen to be conditionally 
positive definite of oz’dei' k=l. 

When the constant r is taken to be zero, each basis function 
is the upper half of a cone, hence the interpolant is not smooth 
at the data points. This particular interpolant is almost the 
"multiconic" method of Duchon [13]. The difference is that the 
multiconic approximation is consistent with the concept of 
conditional positive definiteness of order one. The overall 
approximation takes the form of a linear combination of the basis 
functions plus a constant, with the additional constraint that 
the coefficients satisfy, 

j=i 

This constraint is easily seen to guarantee precision for con- 
stant functions. This special case, along with Theorem 1, is 
convincing evidence that a constant and the corresponding con- 
straint should be included in the MQ approximation. In a practi- 
cal vein, however, my own limited tests have indicated that the 
accuracy of the method is not always helped, and may be hindered, 
by including the constant. 

The "reciprocal" MQ method uses 2F’(t) for the basis 
function, and thus the theorem shows that the coefficient matrix 
for that scheme is positive definite. 
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Theorem 2: Let l=[s/2]-k+2 be a positive integer. Then for any 
function defined on (0,°o ) such that ( -1 ) ^ ^ ( t ) is nonnega- 
tive, nonincreasing, and convex for j=0, 1, ..., 1-2 on (0,cc) (if 
1=1, we require only that it be nonnegative and nonincreasing), 
it follows that F(\/t) is in Pj^(R^). 

Theorem 3: Suppose F’ is completely monotonic but not constant 
on (0, co), F is continuous on C0,«5) and positive on (0,-c). Then 
for any distinct vectors x^, , . . . , x^ in (s arbitrary) 

(-D’^'^detiFC I Ix^-x-^ I | 2 )} > 0 . 

This last theorem proves the con.jecture about the coeffi- 
cient matrix in the MQ method. However, the theorem also yields 
similar results for interpolation by other sets of radial basis 
functions. Some of those methods are known by other results to 
always lead to nonsingular systems of equations, such as the 
"thin plate splines" of Duchon [12-13] (see also Harder and 
Desmarais [21] and Meinguet [37-40]), which are known to be 
nonsingular because of their semi-Hilbert space setting. Others, 
such as the basis function log( 1 -*- )| x-x"^ I !^) suggested by Dyn [16] 
are also seen to be positive definite of order k=l. 

Another result of interest in specific applications has been 
obtained by Hardy and Nelson [26] . This result shows that the MQ 
method has a basis for approximation of geodetic and gravita- 
tional anomalies. The connection comes through the representa- 
tion of the disturbing potential at point i due to anomalous 
disturbances as a three dimensional integral, 

= / d. p dV , 
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where is the distance from the point i, p is the Laplacian of 
the disturbing potential, and V is the volume of interest. What 
Hardy and Nelson showed is that the MQ method can be interpreted 
as a quadrature approximation to the above integral, one that is 
required to yield the correct result at certain (the data) 
points. While this particular result may say little about the 
scheme as a general interpolation scheme, it is interesting that 
such an interpretation exists for one of the early the uses of 
the method . 

The role of the parameter r in the MQ method has long puz- 
zled investigators. The original interpretation given by Hardy 
was a three dimensional one: His applications were actually in 
3-space, and this value simply represented the z coordinate of 
the locations of the disturbing (point) potential. Given that 
the method performs so well, it seems likely that the MQ approxi- 
mation can be described as approximation in some Sobolev-like 
space similar to that for the multiconic method. A recent devel- 
opment by Madych and Nelson [35] gives this result. The MQ 
method does minimize a certain pseudonorm involving a weighted 
integral. Further, the result applies to other similar schemes, 
and thus may also lead to ways of deriving other approximations 
with desired properties through explicit minimization of weighted 
pseudonorms. If practical for computational purposes, such re- 
sults could have far reaching implications in applied scattered 
data approximation. 

2 . 2 Multiquadric Method : Practice 

In addition to the developments on the mathematical aspects 
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of the MQ method, some progress has been made in attempting to 
solve such systems of equations by iterative means. A series of 
investigations by Dyn and coworkers [15-16] have studied the 
problems of fitting scattered data with linear combinations of 
radial basis functions. The basic idea comes from the fact that 
thin plate splines have basis functions which are the fundamental 
solution of the biharmonic equation. The MQ basis functions are 
solutions of Laplace’s equation (in three dimensions). Thus, 
when certain finite difference approximations to the iterated 
Laplacian are applied to the equations, the resulting equations 
tend to have a large diagonal term, which then yields a system of 
equations amenable to solution by iterative schemes. This pro- 
cess may be thought of as applying a conditioning operator to the 
system of equations. 

The idea of the conditioning operator is to transform that 
part of the system of equations involving A into an equivalent 
one which is better conditioned, perhaps even diagonally 
dominant. In addition to transforming A suitably, the operator 
is constructed to annihilate E. This yields a singular system 
CAa = Cf , which has a unique solution, subject to the constraint 
E’a = 0. 

Construction of the conditioning operator involves triangu- 
lation of the convex hull of the data point locations. Finite 
difference operators are then derived which approximate the nec- 
essary derivatives on the basis of the function behavior at the 
vertices of the triangles. Some care is necessary to ensure that 
the operators annihilate E, which corresponds to annihilating all 
polynomials of degree <m on the given set of data points. 
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The results of applying these ideas to only a limited number 
of irregular, but "quasi-regular" grids are reported in [16]. 
Nonetheless, the results are very encouraging, with the condition 
number of the matrix being decreased by factors of up to 200 or 
more for the MQ method with up to 121 data points. In addition 
to the MQ method, the ideas are applicable to approximation by 
thin plate splines and other radial functions , and these are 
reported on as well. 

One further interesting aspect of iterative schemes based on 

this conditioning is that for the MQ method, shifted logarithm 
2 2 

(log(d + r ) ) basis functions, and shifted thin plate spline 

2 2 2 2 

((d + r )log(d + r ) ) basis functions certain nice spectral 

properties of the matrix CA seem to occur. In particular, it was 
noted that computationally, the "rough" eigenvectors correspond 
to the smaller eigenvalues. Certain iterative schemes will re- 
move those components quickly, leaving the surface corresponding 
to the approximate solution (before iteration to convergence) as 
a smooth one. Thus, terminating the iterative scheme prematurely 
corresponds to a smoothing scheme. Numerical evidence and exam- 
ples are given by Dyn, Levin, and Rippa [16]. See [14] in this 
Proceedings for more recent results. 

3.0 Surfaces with Tension or Constraints 

In the univariate case, development of curves with tension 
parameters, constrained approximation, and monotone and convex 
approximation is at a high level. In each situation, several 
reasonable schemes for obtaining such approximations are readily 
available. In two or more independent variables, for scattered 
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data, the situation is not nearly so well developed. Indeed, the 
idea of exactly what characterizes monotone behavior of scattered 
data has not yet been clearly given. Nonetheless, there have 
been some noteworthy developments in the general area of surfaces 
with tension and surfaces satisfying certain constraint relation- 
ships, both from theoretical and computational points of view. 

3 . 1 Tension 

A generalization of splines under tension to the scattered 
data case was given by Nielson and Franke [45] . The basic idea 
is that of approximation of the the surface by a "finite element" 
method, such as proposed by numerous authors, among them Dooley 
[10], Akima [1], Lawson [31], Nielson [43], and LeMehaute [33]. 
The first step of this process consists of triangulation of the 
x-y data in some reasonable way (e.g. , using the max-min angle 
criterion). A certain "finite element" function is assumed over 
each triangle. Ordinarily one wants a smooth surface, so at 
least functions are usually used. Depending on the element 

chosen, certain derivatives must be estimated from the data. It 
is this stage of the process which is crucial to the effective- 
ness of the method (see Nielson and Franke [44]). 

The incorporation of tension into the approximation is 
achieved in the following manner. Let the vertices (x-y data 
pairs) of the triangulation be denoted by the set {V^}, an index 
set of triangle edges joining V. and V. by N , and the edges of 

1 J © 

the triangulation by the set {e. .;ij€N }. The ordering of the 

1 J © 

edges relative to the ordering of the data points is unimportant, 
but each edge must appear only once in the set. Let E be the 
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union of t,he set of edges. Define the set of functions C(E) as 
those which are over the convex hull of the {V^}, restricted 
to E. Let the vector <a. .> be given, with nonnegative compo- 

+ Vi 

nents, each corresponding to a tension parameter for the ij ^ 
edge. Now define the pseudonorm 



01 



F) 



= S 

i.j ^N 



[(52F 2 
’• ’-(9e2 > 
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de. . 

1.3 



1.3 



There is a unique minimiser of S^(F) over all curve networks in 
C(E) which interpolate the data. That minimizer is the function 
which is a piecewise Hermite exponential (a linear combination of 
1, u, exp(a. .u), expl-a. .u), where u represents Euclidean dist- 

-L -1- vJ 

ance along the edge, that takes on prescribed value and deriva- 
tive conditions at the endpoints), and satisfying a certain 
system of (sparse) linear equations for the partial dex’ivatives 
at the vertices. Asymptotic properties are as anticipated. For 
simplicity, take all tension parameters to have the same value, 
and then as tension goes to zero, the curve network approaches 
the curve network minimizing the corresponding functional, as 
developed previously by Nielson [43] . As the tension becomes 
large, the curve network approaches the piecewise linear network 
over the edges. 

After obtaining the curve network over the edges of the 
triangulation, the surface is then completed using a blending 
method on the individual triangles. It is desirable to propogate 
the effects of tension into the interior of the triangles, so the 
previously known cubic blending techniques are inadequate. The 
method used was an extension of the side-vertex scheme (see 
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Nielson [42] ) where the radial projectors were taken to be expo- 



nential Hermite functions. The net effect of the tension becom- 
ing large is that the surface tends toward the the surface which 
is piecewise linear over the triangulation. 

The limiting behavior of the above construction is not the 
physical analog of the limiting behavior of a thin plate under 
tension. In an effort to model the behavior of the thin plat,e 
under tension in two independent variables, Franke [20] developed 
a surface which is the analog of the spline under tension in the 
same way that thin plate splines are the two dimensional analog 
of cubic splines. The "engineering" approach of Harder and 
Desmarais [21] is followed. The idea is to use superposition of 
fundamental solutions for a plate with tension to model the 
displacement under point loads. The fundamental solution (under 
appropriate assumptions about the plate parameters) satisfies 



where, Kq is the modified Bessel function. The solution to the 
interpolation problem is obtained by using this basis function (C 
is taken to be zero) in a global basis function method with m=l. 
Some experimentation was performed including linear polynomial 
terms, since this is more consistent with thin plate splines. 

For tension a = 0, the equation for the thin plate spline is 
attained, and as tension gets large, the equation becomes the 
membrane equation, which has no finite solution under point 



= (5 , 



where a is a tension parameter. The solution is 




dsdt + C 
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loads. Examples are given which demonstrate that as tension is 
increased the surface tends to behave somewhat like a rubber 
sheet under point loads: too thick to be a membrane, but sup- 
porting little displacement away from the data points. 

Construction of surfaces under tension has also been consid- 
ered by Terzopoulous [58], These surfaces are the analog of 
surfaces constructed by Briggs’ method [6], although the paper 
also addresses the solution of the system of equations by multi- 
grid methods. That appears to be quite effective, but will not be 
addressed here. In addition, as with Briggs' scheme, the best 
situation is when the data points all lie on a subset of a 
rectangular grid. The present development has only a under- 
lying surface, which simplifies the calculations for minimization 
of the functional. 

The idea is to minimize the sum of a penalty functional 
measuring closeness of fit to the data and the functional 

P(x,y){r(x,y)(F^^+2F^y+F^y) + ( l-r(x, y ) ) (F^+F^ ) }dA , 

where D is the region of interest in the plane, p(x,y) is the 
"rigidity” of the plate and 0<.r(x,y)<l is the "surface tension". 
A rectangular grid is placed over the region D. While not neces- 
sary in practice, for simplicity it is assumed that the data 
points coincide with grid points. A nonconforming piecewise 
quadratic finite element is assumed, and the equivalent finite 
difference equations for the above functional are derived. These 
equations are solved by multi-grid techniques. Local tension can 
be achieved by allowing the tension parameter to vary over the 
grid points. Discontinuities in derivative and value are also 
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considered. As with Briggs' method, the final form of the ap- 
proximation depends on the grid at which it is evaluated. Unlike 
Briggs’ scheme, which is based on finite difference approxima- 
tions to the plate equation, the underlying surface (although C^) 
would allow evaluation at other than grid points. 

3 . 2 Constraints 

The construction of approximations satisfying inequality 
constraints as well as interpolation conditions has been investi- 
gated by Villalobos [62], and Dubrule and Kostov [11]. The 
former considered a generalization of Laplacian smoothing splines 
(see Section 4.1), while the latter considers only interpolating 
functions. The results are similar, and the latter will be 
described here. The discussion centers on approximation by thin 
plate splines. The functional minimized is the usual thin plate 
functional, but under the given constraints. As had been pre- 
viously shown by others, the solution involves adding basis 
functions at constraint points that are found to be active, and 
is solved as a quadratic programming problem. This method ap- 
pears to be easy and effective. The practical aspects of the 
method are discussed in the companion paper [30], where the 
process is applied using Kriging, which formally includes thin 
plate splines as a special case. Journel [28] discusses the 
incorporation of this and other "soft" information into the 
approximation by Kriging. 

A more general problem and its elegant solution is consid- 
ered by Utreras [61] . Here the problem is that of minimizing the 
thin plate functional, 
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/ (F^ +2F^ +F^ ) dA 

I XX xy yy ' 

subject to interpolation conditions, and positivity conditions on 
a certain subset of the plane, say 
F(x,y) > 0 for (x,y)€ KCR2. 

Here the region K is assumed compact and convex. The existence 
and uniqueness of the solution is proven using elementary means. 
The characterisation of the solution of the problem is in terms 
of the points where the constraints are active. Distribution 
theory and the relationship of the functional minimised to the 
biharmonic equation are used to show that the solution involves 
the usual terms of the thin plate spline approximation plus a 
certain term which serves to enforce the positivity of the func- 
tion over the compact region K. 

We digress to introduce some notation. Suppose that scat- 
tered data is given with f^>0. Let F(x,y) be the solution of the 
constrained problem which interpolates this data. Define the set 
K’ = {(x,y)€ K:F(x,y) = 0}. K’ is compact. Let B^(x,y) = 
d^log d-^, where d^ = + , and d^ = x^+y^ . The 

following theorem, where * denotes convolution, is then proved. 

Theorem; There exist constants and a positive measure fi with 

support contained in K’ such that 

N 

F(x,y) = ^^a^ Bj^(x,y) + lft*BQ(x,y) + Pj^(x,y) , 
k=l 

where p^ is a polynomial of degree <1, and 

N 

k=l 

for any polynomial q of degree <1. 
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Convergence is investigated, and an algorithm is given for 
computing an approximation to the positive thin plate spline. 
The algorithm involves finding an approximation of the set K’ and 
the measure A such that the constraint is satisfied to within 
some tolerance. The approximation is by points with atomic 
measure, and since convolution with atomic measures give point 
evaluation, the solution is approximated by functions of the same 
type as in the Dubrule/Kostov program. The crucial difference is 
that here the constraint set K’ must be completely discovered as 
part of the process, rather than being part of a finite subset 
specified in advance. 

4.0 Smoothing, Least Squares, and Subset Selection 

In many applications the data is obtained by measurement, 
often not to enough accuracy to warrant interpolation. In other 
applications, such as oceanography and remote sensing, the amount 
of data available, even though subject to errors, is much greater 
than is necessary to define the desired surface to the required 
accuracy. In these cases is is necessary to apply some smoothing 
process or to perform least squares approximation with some 
function having far fewer parameters than the number of data 
points. The theory of smoothing splines in several variables is 
well developed. For purposes of contrasting that idea with those 
of least squares and subset selection, a brief discussion of 
smoothing splines will be given. 

4.1 Laplaclau Smoothing Splines 

Laplacian smoothing splines are a generalization of univa- 
riate smoothing splines to several variables along the same 
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direction as thin plate splines for interpolation. They were 

mentioned by Harder and Desmarais [21] in their seminal paper for 

the bivariate case under the physical interpretation of having 

forces applied at the data points through springs with various 

spring constants. This caused the surface to tend toward the 

data points, but the not necessarily to pass through them. 

The mathematical development of Laplacian smoothing splines 

in the general case is given in Wahba and Wendelberger [64] . The 

notation of Section 2.1 for points in s-dimensional space will be 

followed in this discussion. The functional minimized in the 

case of smoothing splines is 
N 

N‘^^[F(x'^)-f^]2/a2 + XJ^(F) 

j = l 

where is the pseudonorm associated with interpolating splines, 

and X is a smoothing parameter which governs the fidelity with 
which the surface fits the given data. Here it has been assumed 
that the errors are uncorrelated and have standard deviation <7^ 
at the point . The case of correlated errors is addressed 
briefly in Wendelberger ’ s thesis [65]. 

Let ct - , 02, , a^) be a multi-index for partial 

differentiation (denoted by F ) , with each integer ^0 and |0t | 

= ••• ■'■‘^s' 

|a|=m-^R® ’ 

In this functional one can think of m as a smoothness parameter, 
the approximating functions being smoother (having higher order 
derivatives) the larger m is. In order for the required func- 
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tions "to be in the appropriate spaces, 
m>s/2 . 



it is necessary that 



The solution of the problem is of the same form as thin 

plate splines, 

N 

S-(x) = t <3„.i(x) , 

J=1 



where the coefficients a^ and the coefficients bg, of ^ 

polynomial of degree m-1, satisfy the system of equations 



N 

j=l ia|<m 



N 

= 0, lal<m. 

i=l 



In the above equations, (5^^ is the Kronecker delta. 

, , ,a oJi ag 

and ( X ) = X. Xr, . . . x 

1 ^ 3 

The smoothing parameter X must be specified before the 
smoothing spline can be computed, and Wahba and Wendelberger 
suggest the use of Generalized Cross Validation (GCV) to select 
the value. Although there is some evidence that GCV leads to 
undersmoothing (see Seaman and Hutchinson [56]), Utreras [60] has 
shown that GCV yields convergence under reasonable conditions. 
GCV can also be used to decide on the smoothness parameter, m, to 
be used, as well. A limitation of the scheme is that it is 
difficult to compute when there are more than 200-300 data 
points, since the problem involves the solution of a system of 
more than N linear equations, although techniques such as those 
of Dyn, described in the previous section, could be useful here. 
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4 . 2 Least Squares Approximation 

Least squares approximation to scattered data is an approxi- 
mation with (presumably many) fewer basis functions than there 
are data points. The problem at hand here, then, is to select 
the basis functions. Previous work has been done using tensor 
product cubic splines, and a number of authors (see [T] , [8], 
[9], for example) have considered the problem. Several computer 
programs are available, and such methods are probably desirable 
for cases where data is somewhat uniformly distributed. In cases 
where the data is of greatly varying density, the use of tensor 
products results in knot locations on a grid, and this may not 
reflect the actual disposition of data points. In fact, there 
could be knots with no data nearby. While such problems are not 
insurmountable, they lead to nonuniqueness of the solution, and 
the minimum norm solution tends to not be aesthetically pleasing. 

For varying density of data points it seems desirable to 
have flexibility in knot placement, and this leads to the idea of 
least squares approximation by thin plate splines. The MQ ap- 
proximations also could be used, however the discussion will be 
in terms of thin plate splines, and the points at which these 
basis functions are centered will be referred to as "knot" 
points. Treatment of the knot point locations as parameters in 
the minimization process is possible, and has been reported upon 
by Schmidt [53]. The paper is brief, with few details of the 
algorithm being given. The initial knot configuration was taken 
to be of tensor product form, which may be apparent from the 
final configuration of knots in the examples given. The overall 
minimization process is a large nonlinear one, and if the 



20 



problems of one dimension carry over as one might expect, may be 
complicated in that knots may coalesce and the solution may not 
be unique. In addition it is likely true that the objective 
function may have many local minima, so an algorithm to search 
for a better local minimum, or to avoid a poor local minimum 
would be necessary. 

The problem with treating the knot locations as parameters 
presently seems to be intractable mathematically, and for many 
knots is computationally expensive, with results obtained being 
of questionable quality. A different point of view on the prob- 
lem is considered in McMahon [36] , and a summary of his approach 
and results will be given. 

The main problem in the process is the selection of knots, 
for once these have been decided upon, what remains is to solve 
the equations in the least squares sense, in principle an unin- 
teresting task. If the selection of knot locations is to be 
decoupled from the least squares process, some assumption must be 
made order to have an algorithm for selection of the knots. The 
assumption in this case is that the independent variable data 
reflects something about the behavior of the dependent variable. 
For example, perhaps the density of data point locations is 
dependent on the curvature of the surface, or more broadly, if 
the function is changing behavior rapidly the density of data 
points is great, whereas low density indicates slowly changing 
behavior. This assumption is not universally satisfied in prac- 
tice, for some data is taken based on accessibility (along roads 
in rugged areas, for example) or other nonbehavioral criteria. 

The assumption of data density indicating local behavior of 
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the surface has a companion assumption that each data point is in 
some sense equally important in defining the function, regardless 
of whether it has a nearby neighbor or not. This leads to the 
idea of "equal representation" for each data point by a knot 
point. This means that each data point should be "close" to a 
knot point, and that each knot point should "represent" about the 
same number of data points. These two ideas are crucial to the 
knot selection algorithm. First, it is desirable to minimize the 
sum of the distances squared from each data point to the nearest 
knot point. Let the knot point locations be given by 

0 J 

j=l, 2, . . . , K. Then the function to be minimized is 

N 






k=l 



This process has a default Dirichlet tesselation with respect to 
the knot points, with each data point "belonging" to some knot 
point by virtue of the Dirichlet tile to which it belongs. When 
a data point lies on a tile boundary, some determination of which 
it belongs to, or whether it is shared must be assumed. Local 
minima for the above function are easily characterized; At a 
local minimum, each knot point (X.,J?.) is the centroid of the 
data points in its tile. This leads to a nice algorithm for 
iteration to a local minimum, by repeated computation of the 
centroid of the data points in the Dirichlet tiles. 

This algorithm works very well for obtaining a local minimum 
2 

of the function GN . Unfortunately the function is rife with 
local minima, so that finding a desirable one depends on the 
proper initial guess. Here the second idea regarding "equal 
representation" comes into play. Since each data point is 
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assumed equally important, it is reasonable to expect that the 

Dirichlet tile for each knot should contain about the same number 

of data points. This leads to a heuristic for a new initial 

guess at knot locations once a local minimum for GN has been 

found. It is still considered desirable to be at a local minimum 

of GN , but the idea of finding a global minimum of GN is 

abandoned for the "equal representation" idea. Once a local 

minimum has been found, a new measure of "goodness" is computed. 

t h 

Let N ^ be the number of points in the Dirichlet tile for the 

knot point ( ii . , ^ . Let 

J 0 

K 

D = ^(Nj-N/K)2 , 
j = l 

a measure of the "equal representation" for a particular knot 

configuration. It was then attempted to determine knot locations 

2 

which achieve a local minimum of GN and which also yield a small 

value, hopefully a minimum, of D. 

The general idea of the algorithm is to "nudge” the knot 

2 

locations from a local minimum of GN toward a configuration with 
smaller D value. The rationale for this is to attempt to move 
the tile boundaries across data points, yielding a more equitable 
distribution. The search for the minimum of D could turn out to 
be rather extensive, and for a large number of points with a 
moderately large number of knot points, the computational effort 
can be excessive. End results are still somewhat dependent on 
the initial guess, although the results generally tend to look 
quite reasonable. The program incorporates the option of inter- 
nally generated ( quasi-gridded ) or user input initial guesses. 
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Examples are given which illustrate rather well the ability 
of the scheme to select knot locations which reflect the under- 
lying density of the data. Actual surface fitting and comparison 
with two other methods, the Laplacian smoothing splines of Wahba 
and Wendelberger , and the tensor product bicubic Hermite method 
due to Foley [18], are reported upon. One example illustrates 
the failure of the fitting scheme to accurately model the surface 
when the data density does not reflect the underlying behavior of 
the surface, as was assumed. 

4 . 3 Subset Selection 

To my knowledge, the first investigation into the use of a 
subset of scattered data points upon which to base an interpolant 
to be used as an approximation for the entire data set was by 
Pickrell [47] . The application guided many of the ideas involved 
in the thesis and they are not universally applicable. 
Pickrell' s goal was to model underwater terrain from a large 
number of measured depths, which were quite accurate, the most 
important being those which were on ridges or other shallow 
areas, since the resulting approximation could be used to gener- 
ate charts for navigation. 

The idea was to select as small a subset of points as pos- 
sible such that the interpolant for these points yielded a sur- 
face that satisfied a certain error tolerance at all the other 
data points. An iterative scheme was used. Beginning with an 
initial guess, either taken "uniformly" spaced over the region of 
interest, or by perusal of the data for "critical" points such as 
high and low values or areas of sharp gradients. Call the selec- 
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ted points the "model" points. The interpolant (the MQ method 
was used; parameter r=0 generally performed best) was constructed 
and the deviation at other data points computed. Based on this 
information, other points at which the deviation was large, or 
one from a group of nearby points where deviations of the same 
sign occurred, would be included in the set of model points. 
Points with small coefficients might be eliminated from tne set 
of model points, as well. The investigation did not yield an 
algorithm suitable for approximation with no intervention by the 
user, since this was performed interactively by Pickrell. In the 
examples reported, the required number of model points was gen- 
erally on the order of 10-20% of the total number of points. 
However, only limited consideration was given to the handling of 
very large data sets which would require local application of the 
ideas with a scheme for .joining the pieces together. 

These deficiencies and other matters were addressed by 
Schiro and Williams [52]. In addition to automating the model 
point selection process, he subdivided the region into kernel 
groups based on a measure of homogeneity of the data. As a 
starting point the region of interest was divided into a rectang- 
ular grid of cells of size equal to the minimum to be considered. 
The mean and standard deviation of the function values (again, 
ocean depth data) were computed for each cell. Then, the cells 
were considered in turn as a base cell, and contiguous cells 
(with larger coordinates) having mean values within one standard 
deviation of the mean for the base cell were combined to form a 
larger rectangular cell, when possible. This larger group of 
cells is called a kernel group. Thus, each kernel group in the 
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final subdivision was made up of one or more contiguous cells in 
the original subdivision, all of which have similar mean behavior 

The selection of model points for each kernel group is done 
without user intervention. The initial set is chosen based on 
deviation from the mean value. (Of course, if all points are 
within the tolerance of the mean value, the mean value is used as 
the approximation for the kernel group. ) The idea is similar to 
Pickrell’s: points with large deviation are selected to be added 
to the set of model points, with the proviso that points closer 
than a certain distance cannot be selected on the same iteration. 
This avoids the problem of adding several points very close 
together on the same iteration. The process is continued until 
the deviations are below a specified tolerance. 

The overall surface is made continuous by blending adjacent 
surfaces when within a certain distance from boundaries of the 
kernel groups. This is done via Hermite cubics to obtain a 
smooth transition. The MQ method with r=0 was used to fit the 
differences, data minus mean value, and this was observed to have 
a beneficial effect in terms of the magnitude of the coefficients 
in the representation. 

Another approach to subset selection was taken by Bozzini , 
deTisi, and Lenarduzzi [5]. Here an attempt was made to deter- 
mine a subset of points to be used to compute the approximation 
by doing local computations to determine whether a particular 
point has a significant influence on defining the surface. A 
description of the ideas will be given. The first step is to 
partition the region of interest into (probably overlapping) 
regions of given shape (e.g., circular disks) by taking the 
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region as large as possible so that 



— \l lf(x)-f . I dx <. K. , for some constant K. . 

\ xV ^ ^ ^ ^ ^ 



The value of specifies a kind of homogeneity of f(x) over R^. 
Then, for a weighting function <i>^(x), let 



3. 

1 



]f(x)-f.| 4(x) dx // i>(x) dx 

«i 



This value serves as measure of the level of homogeneity in the 
region , weighted by . Then a measure of the behavior of f(x) 
at point i relative to point j (point j also in ) is given by 



1 J 



I f _f I 

■^i 



/<Ji 






These values are averaged over the points in to obtain s^. 
The value of 3^ is a measure of the what the authors call the 
"importance" of the point i in the data set, and hence to the 
definition of the approximating surface. The integrals in the 
definition are approximated with the obvious quadratures in the 
application. The point with the largest corresponding value of 
is chosen. The process can then be repeated from the computa- 
tion of the to obtain more points, or one can choose the 
subset with the largest values from the initial calculation. 
According to the authors, the method is not to sensitive to 
deletion of a point from the set. Some examples are given to 
illustrate approximation of surfaces from subsets of a given set 
of points. 
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